# Clear memory
rm(list = ls())

# Load packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  maps, data.table, janitor, broom, haven, sf,
  statar, fixest, patchwork, tidyverse, tidylog, showtext, rmapshaper, 
  maps
)
options("tidylog.display" = NULL)

font_add_google("Lato")
showtext_auto()

#################################
#################################
## LOAD DATA
#################################
#################################

state_fips <- state.fips %>%
  distinct(fips, .keep_all = TRUE)

dup_names <- str_sub(state_fips$polyname, 1, str_locate(state_fips$polyname, ":")[, 1] - 1)

state_fips$polyname[!is.na(dup_names)] <- dup_names[!is.na(dup_names)]

state_sdf <- st_as_sf(maps::map("state", plot = FALSE, fill = TRUE)) %>%
  rename(polyname = ID) %>%
  full_join(state_fips) %>%
  st_transform(5070)
state_sdf_simpler = rmapshaper::ms_simplify(state_sdf, keep_shapes = TRUE, keep = .05)

county_fips <- county.fips
counties_sdf <- st_as_sf(maps::map("county", plot = FALSE, fill = TRUE)) %>%
  rename(polyname = ID) %>%
  mutate(polyname = ifelse(polyname == "south dakota,oglala dakota", "south dakota,shannon", polyname)) %>%
  full_join(county_fips) %>%
  mutate(fips = ifelse(fips == 46113, 46102, fips)) %>%
  mutate(fips = ifelse(polyname == "florida,okaloosa", 12091, fips)) %>%
  mutate(fips = ifelse(polyname == "virginia,accomack", 51001, fips)) %>%
  mutate(fips = ifelse(polyname == "texas,galveston", 48167, fips)) %>%
  mutate(fips = ifelse(polyname == "louisiana,st martin", 22099, fips)) %>%
  mutate(fips = ifelse(polyname == "north carolina,currituck", 37053, fips)) %>%
  mutate(fips = ifelse(polyname == "washington,pierce", 53053, fips)) %>%
  mutate(fips = ifelse(polyname == "washington,san juan", 53055, fips)) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  st_transform(5070)
counties_sdf_simpler = rmapshaper::ms_simplify(counties_sdf, keep_shapes = TRUE, keep = .05)

non_attainment <- fread("data/productivity_panel.csv") %>%
  filter(year >= 1991 & year <= 1997) %>%
  group_by(fips, industry) %>%
  mutate(total_nonattainment = sum(nonattainment, na.rm = T)) %>%
  ungroup() %>%
  mutate(nonattainment_1991 = nonattainment > 0 & year == 1991) %>%
  mutate(nonattainment_1997 = nonattainment > 0 & year == 1997) %>%
  group_by(fips) %>%
  mutate(nonattainment_1991 = max(nonattainment_1991)) %>%
  mutate(ever_nonattainment = max(nonattainment_1997) > 0) %>%
  ungroup() %>%
  filter(year == 1997) %>%
  mutate(fips = str_pad(fips, 5, "left", "0")) %>%
  select(fips, year, industry, nonattainment_1991, nonattainment_1997, ever_nonattainment, emp) %>%
  filter(industry == "Manufacturing")

payroll <- fread("data/productivity_panel.csv") %>%
  group_by(fips, year) %>%
  mutate(emp = emp / sum(emp)) %>%
  filter(industry == "Manufacturing") %>%
  mutate(fips = str_pad(fips, 5, "left", "0")) %>%
  mutate(pay = pay / emp / 1000) %>%
  filter(year == 1997 | year == 1991) %>%
  mutate(year2 = year) %>%
  pivot_wider(
    names_from = "year",
    names_prefix = "emp_",
    values_from = "emp"
  ) %>%
  group_by(fips) %>%
  mutate(emp_1991 = lag(emp_1991)) %>%
  filter(year2 == 1997) %>%
  rename(year = year2) %>%
  ungroup() %>%
  arrange(desc(emp_1991)) %>%
  mutate(rank_1991 = row_number()) %>%
  arrange(desc(emp_1997)) %>%
  mutate(rank_1997 = row_number()) %>%
  mutate(emp_change = (emp_1997 - emp_1991) / emp_1991) %>%
  arrange(fips) %>%
  filter(industry == "Manufacturing")

plot_sdf <- counties_sdf_simpler %>%
  left_join(non_attainment) %>%
  left_join(payroll, by = c("fips", "year", "industry")) %>%
  st_transform(5070)

plot_sdf$post_nonattainment[is.na(plot_sdf$post_nonattainment)] <- FALSE
plot_sdf$nonattainment_1991[is.na(plot_sdf$nonattainment_1991)] <- FALSE
plot_sdf$ever_nonattainment[is.na(plot_sdf$ever_nonattainment)] <- FALSE

#################################
#################################
## NONATTAINMENT MAPs
#################################
#################################

ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = ever_nonattainment),
    col = "#DDDDDD",
    size = 0.2
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.3,
    fill = "transparent"
  ) +
  scale_fill_manual(
    name = "",
    labels = c("In Attainment in 1997", "In Nonattainment in 1997"),
    values = c("#ffffef", "#b2182b")
  ) +
  theme_minimal() +
  theme(
    axis.title = element_text(size = 1),
    panel.grid = element_blank(),
    title = element_text(size = 9),
    axis.text = element_text(size = 1, colour = "#FFFFFF"),
    legend.position = "bottom",
    legend.title = element_text(size = 16),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1, "cm"),
    legend.text = element_text(size = 24)
  )


ggsave(
  "output/nonattainment_ever_map.eps",
  width = 20,
  height = 12,
  dpi = 175
)
ggsave(
  "output/nonattainment_ever_map.pdf",
  width = 20,
  height = 12
)


# Black and white for publication 
ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = ever_nonattainment),
    col = "#DDDDDD",
    size = 0.2
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.3,
    fill = "transparent"
  ) +
  scale_fill_manual(
    name = "",
    labels = c("In Attainment in 1997", "In Nonattainment in 1997"),
    values = c("white", "black")
  ) +
  theme_minimal() +
  theme(
    axis.title = element_text(size = 1),
    panel.grid = element_blank(),
    title = element_text(size = 9),
    axis.text = element_text(size = 1, colour = "#FFFFFF"),
    legend.position = "bottom",
    legend.title = element_text(size = 16),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1, "cm"),
    legend.text = element_text(size = 24)
  )

ggsave(
  "output/figure-2-nonattainment_ever_map_bw.eps",
  width = 20,
  height = 12,
  dpi = 175
)




#################################
#################################
## AP3 EXAMPLE
#################################
#################################

fips = fread("raw/modified-ap3/fips.csv") 

fips_order = fread("raw/modified-ap3/fips.csv") 
fips_order = fips_order %>% 
  mutate(order = seq(1:nrow(fips_order))) %>%
  rename(fips = V1)

ap3_sdf <- fread("data/modified-ap3/NOX_transport_L_1997_prime_aged.csv") %>%
  as_tibble()  %>%
  mutate(order = seq(1:nrow(fips_order))) %>%
  left_join(fips_order, by = c("order")) %>%
  select(fips, starts_with("V")) %>%
  rename(fips_origin = fips) %>%
  pivot_longer(cols = V1:V3109, names_to = "name", values_to = "conc") %>%
  mutate(order = rep(seq(1:nrow(fips_order)), times = 3109)) %>%
  left_join(fips_order, by = c("order")) %>%
  rename(fips_destination = fips) %>%
  rename(fips = fips_origin) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips))

la_ap3_sf <- ap3_sdf %>%
  filter(fips == "06037") %>%
  select(-fips) %>%
  rename(fips=fips_destination) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0"))  %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  mutate(main = fips == "06037") %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips))


stl_ap3_sf <- ap3_sdf %>%
  filter(fips == "29189") %>%
  select(-fips) %>%
  rename(fips=fips_destination) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0"))  %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  mutate(main = fips == "29189") %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips))


la_ap3_sf <- counties_sdf_simpler %>%
  left_join(la_ap3_sf) %>%
  st_transform(5070)

stl_ap3_sf <- counties_sdf_simpler %>%
  left_join(stl_ap3_sf) %>%
  st_transform(5070)

ggplot() +
  geom_sf(
    data = la_ap3_sf %>% drop_na(conc),
    aes(fill = conc * 1e6 + 5e-4, color = main, size = main)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#888888", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradient(
    name = "Change in\nParticulate\nConcentrations",
    low = "#ffffff", high = "#b2182b",
    space = "Lab", na.value = "#555555",
    aesthetics = "fill",
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#111111",
      frame.colour = "#111111",
      frame.linewidth = 2,
      labels = waiver()
    ),
    trans = "log10",
    breaks = c(1e-3, 1e-1, 10),
    labels = c("0.001", "0.1", "10")
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.97, 0.30)
  )

ggsave(
  "output/los_angeles_concentrations.eps",
  width = 20,
  height = 10,
  dpi = 175
)
ggsave(
  "output/los_angeles_concentrations.pdf",
  width = 20,
  height = 10
)



ggplot() +
  geom_sf(
    data = stl_ap3_sf %>% drop_na(conc),
    aes(fill = conc * 1e6 + 5e-4, color = main, size = main)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#888888", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradient(
    name = "Change in\nParticulate\nConcentrations",
    low = "#ffffff", high = "#b2182b",
    space = "Lab", na.value = "#555555",
    aesthetics = "fill",
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#111111",
      frame.colour = "#111111",
      frame.linewidth = 2,
      labels = waiver()
    ),
    trans = "log10",
    breaks = c(1e-3, 1e-1, 10),
    labels = c("0.001", "0.1", "10")
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.97, 0.30)
  )

ggsave(
  "output/st_louis_concentrations.eps",
  width = 20,
  height = 10,
  dpi = 175
)

ggsave(
  "output/figure-1-a.eps",
  width = 20,
  height = 10,
  dpi = 175
)

ggsave(
  "output/st_louis_concentrations.pdf",
  width = 20,
  height = 10
)

fips = fread("raw/modified-ap3/fips.csv") 

fips_order = fread("raw/modified-ap3/fips.csv") 
fips_order = fips_order %>% 
  mutate(order = seq(1:nrow(fips_order))) %>%
  rename(fips = V1)


ap3_sdf <- fread("data/modified-ap3/NOX_md_L_1997_prime_aged.csv") %>%
  as_tibble() %>%
  mutate(order = seq(1:nrow(fips_order))) %>%
  left_join(fips_order, by = c("order")) %>%
  select(fips, starts_with("V")) %>%
  rename(fips_origin = fips) %>%
  pivot_longer(cols = V1:V3109, names_to = "name", values_to = "conc") %>%
  mutate(order = rep(seq(1:nrow(fips_order)), times = 3109)) %>%
  left_join(fips_order, by = c("order")) %>%
  rename(fips_destination = fips) %>%
  rename(fips = fips_origin) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips))

la_ap3_sf <- ap3_sdf %>%
  filter(fips == "06037") %>%
  select(-fips) %>%
  rename(fips=fips_destination) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0"))  %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips))  %>%
  mutate(main = fips == "06037") %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips))


stl_ap3_sf <- ap3_sdf %>%
  filter(fips == "29189") %>%
  select(-fips) %>%
  rename(fips=fips_destination) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0"))  %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  mutate(main = fips == "29189") %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips))


la_ap3_sf <- counties_sdf_simpler %>%
  mutate(fips = ifelse(fips == "46102", "46113", fips)) %>%
  left_join(la_ap3_sf) %>%
  st_transform(5070) %>%
  mutate(main = fips == "06037" | fips == "29189")

stl_ap3_sf <- counties_sdf_simpler %>%
  mutate(fips = ifelse(fips == "46102", "46113", fips)) %>%
  left_join(stl_ap3_sf) %>%
  st_transform(5070)

# difference in damages
la_ap3_sf$conc <- la_ap3_sf$conc - stl_ap3_sf$conc


ggplot() +
  geom_sf(
    data = la_ap3_sf %>% drop_na(conc),
    aes(fill = conc * 1e3, color = main, size = main)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#888888", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = expression(paste("Change in\nDamages\nper Capita ($)")),
    colors = rev(c("#b2182b", "#f4a582", "#ffffff", "#4393c3", "#2166ac")),
    na.value = "#888888",
    aesthetics = "fill",
    limits = c(-.04, .04),
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    values = scales::rescale(c(-0.04, -0.01, 0, 0.01, .04)),
    oob = scales::squish,
    breaks = scales::pretty_breaks()
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.97, 0.30)
  )

ggsave(
  "output/stl_to_la_damages.eps",
  width = 20,
  height = 10,
  dpi = 200
)

ggsave(
  "output/figure-1-b.eps",
  width = 20,
  height = 10,
  dpi = 200
)
ggsave(
  "output/stl_to_la_damages.pdf",
  width = 20,
  height = 10
)

# ========================
# Share remaining in nonattain
# ========================

fread("data/productivity_panel_time_series.csv") |>
  filter(treated == 1 & year <= 2001) |>
  group_by(year) |>
  summarise(share_still_nonattain = mean(sum_nonattainment) * 100) |>
  filter(share_still_nonattain > 0.5) |>
  ggplot() +
  geom_line(aes(x = year, y = share_still_nonattain), size = 1.5, color = "black") +
  geom_point(aes(x = year, y = share_still_nonattain), size = 6, color = "black") +
  scale_shape(solid = T) +
  theme_minimal() +
  theme(
    legend.position = "none",
    title = element_text(size = 24),
    axis.text.x = element_text(size = 30), axis.text.y = element_text(size = 30),
    axis.title.x = element_text(size = 30), axis.title.y = element_text(size = 30),
    panel.grid.minor.x = element_blank(), panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(), panel.grid.major.x = element_blank(),
    axis.line = element_line(colour = "black"), axis.ticks = element_line(),
    axis.title.y.right = element_text(vjust = 2), axis.title.y.left = element_text(vjust = 2),
    axis.title.x.top = element_text(vjust = 1), plot.margin = margin(15, 15, 15, 15)
  ) +
  labs(
    y = "Share of Counties\nStill in Nonattainment",
    x = ""
  ) +
  coord_cartesian(ylim = c(0, 100))

ggsave("output/nonattainment_remaining.pdf", width = 8, height = 8)
ggsave("output/figure-e1-nonattainment_remaining.eps", width = 8, height = 8)

############################

# first-best emission prices map nh3

# employment in the model is normalized to 1 so need to multiply by total employment because we are summing over all people harmed
employment <- fread("data/payroll_cf.csv") |>
  pull(employment) |>
  sum()

# need 1000*employment adjustment because model wages is in thousands
# and the price is back computed using: ( w .* L ./ e ) .* (ξ ./ (γ .* (1 .- sum(ξ))))'
# w is in thousands, and L is normalized to one so the price units are not correct for plotting (although internally consistent in the model)
# so we need to multiply by 1000 and the actual employment to get the dollar version of w * L
# and thus the dollars per ton for the emission price
adjustment <- employment * 1000

fips <- fread("data/simulation_fips_crosswalk.csv")

non_attainment <- fread("data/productivity_panel.csv") %>%
  filter(year >= 1986 & year <= 1997) %>%
  group_by(fips, industry) %>%
  mutate(total_nonattainment = sum(nonattainment, na.rm = T)) %>%
  ungroup() %>%
  mutate(nonattainment_1991 = nonattainment > 0 & year == 1991) %>%
  mutate(nonattainment_1997 = nonattainment > 0 & year == 1997) %>%
  group_by(fips) %>%
  mutate(nonattainment_1991 = max(nonattainment_1991)) %>%
  mutate(ever_nonattainment = max(nonattainment_1997) > 0) %>%
  ungroup() %>%
  filter(year == 1997) %>%
  mutate(fips = str_pad(fips, 5, "left", "0")) %>%
  distinct(fips, nonattainment_1997)

# Load the county data from the maps package
state_fips <- state.fips %>%
  distinct(fips, .keep_all = TRUE)

dup_names <- str_sub(state_fips$polyname, 1, str_locate(state_fips$polyname, ":")[, 1] - 1)

state_fips$polyname[!is.na(dup_names)] <- dup_names[!is.na(dup_names)]


simulation_results <- read.csv("data/counterfactual/welfare_opt-opt_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(nh3_price_change = cf_emissions_price_nh3 / base_emissions_price_nh3) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  filter(industry == "manufacturing") |>
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment) |>
  mutate(
    cf_emissions_price_nh3 = cf_emissions_price_nh3 * adjustment,
    base_emissions_price_nh3 = base_emissions_price_nh3 * adjustment
  ) |>
  select(fips, nh3_price_change, nonattainment_1997, first_best_emissions_price_nh3 = cf_emissions_price_nh3) %>%
  mutate(first_best_emissions_price_nh3 = ifelse(is.infinite(first_best_emissions_price_nh3), 0, first_best_emissions_price_nh3))

simulation_results_base <- read.csv("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  filter(industry == "manufacturing") |>
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment) |>
  mutate(
    cf_emissions_price_nh3 = cf_emissions_price_nh3 * adjustment,
    base_emissions_price_nh3 = base_emissions_price_nh3 * adjustment
  ) |>
  select(fips, nonattainment_1997, no_reg_emissions_price_nh3 = cf_emissions_price_nh3, nonattain_emissions_price_nh3 = base_emissions_price_nh3) %>%
  mutate(no_reg_emissions_price_nh3 = ifelse(is.infinite(no_reg_emissions_price_nh3), 0, no_reg_emissions_price_nh3)) %>%
  mutate(nonattain_emissions_price_nh3 = ifelse(is.infinite(nonattain_emissions_price_nh3), 0, nonattain_emissions_price_nh3))

simulation_results <- simulation_results |>
  left_join(simulation_results_base) |>
  mutate(
    first_best_emissions_price_nh3 = first_best_emissions_price_nh3 - no_reg_emissions_price_nh3,
    nonattain_emissions_price_nh3 = nonattain_emissions_price_nh3 - no_reg_emissions_price_nh3,
  )

plot_sdf <- counties_sdf_simpler %>%
  left_join(simulation_results) %>%
  replace_na(list(nonattainment_1997 = FALSE)) %>%
  st_transform(5070) |>
  st_simplify(dTolerance = 1e3)

ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = (first_best_emissions_price_nh3 + 1) / 1000, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "NH3\nFirst-Best\nEmissions\nPrice\n(thousand $)",
    colors = rev(c("#b2182b", "#f4a582", "#fddbc7", "#ffffff")),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    trans = "log10",
    labels = scales::label_dollar(),
    limits = c(NA, 100),
    breaks = c(.01, .1, 1, 10, 100)
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.96, 0.30)
  )

ggsave(
  paste0("output/first_best_prices_nh3.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)
ggsave(
  paste0("output/figure-e5-a-first_best_prices_nh3.eps"),
  width = 20,
  height = 10,
  dpi = 175
)

# first-best emission prices map so2

# need 1000*employment adjustment because model wages is in thousands
# and the price is back computed using: ( w .* L ./ e ) .* (ξ ./ (γ .* (1 .- sum(ξ))))'
# w is in thousands, and L is normalized to one so the price units are not correct for plotting (although internally consistent in the model)
# so we need to multiply by 1000 and the actual employment to get the dollar version of w * L
# and thus the dollars per ton for the emission price

non_attainment <- fread("data/productivity_panel.csv") %>%
  filter(year >= 1986 & year <= 1997) %>%
  group_by(fips, industry) %>%
  mutate(total_nonattainment = sum(nonattainment, na.rm = T)) %>%
  ungroup() %>%
  mutate(nonattainment_1991 = nonattainment > 0 & year == 1991) %>%
  mutate(nonattainment_1997 = nonattainment > 0 & year == 1997) %>%
  group_by(fips) %>%
  mutate(nonattainment_1991 = max(nonattainment_1991)) %>%
  mutate(ever_nonattainment = max(nonattainment_1997) > 0) %>%
  ungroup() %>%
  filter(year == 1997) %>%
  mutate(fips = str_pad(fips, 5, "left", "0")) %>%
  distinct(fips, nonattainment_1997)


simulation_results <- read.csv("data/counterfactual/welfare_opt-opt_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(so2_price_change = cf_emissions_price_so2 / base_emissions_price_so2) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  filter(industry == "manufacturing") |>
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment) |>
  mutate(
    cf_emissions_price_so2 = cf_emissions_price_so2 * adjustment,
    base_emissions_price_so2 = base_emissions_price_so2 * adjustment
  ) |>
  select(fips, so2_price_change, nonattainment_1997, first_best_emissions_price_so2 = cf_emissions_price_so2) %>%
  mutate(first_best_emissions_price_so2 = ifelse(is.infinite(first_best_emissions_price_so2), 0, first_best_emissions_price_so2))

simulation_results_base <- read.csv("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  filter(industry == "manufacturing") |>
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment) |>
  mutate(
    cf_emissions_price_so2 = cf_emissions_price_so2 * adjustment,
    base_emissions_price_so2 = base_emissions_price_so2 * adjustment
  ) |>
  select(fips, nonattainment_1997, no_reg_emissions_price_so2 = cf_emissions_price_so2, nonattain_emissions_price_so2 = base_emissions_price_so2) %>%
  mutate(no_reg_emissions_price_so2 = ifelse(is.infinite(no_reg_emissions_price_so2), 0, no_reg_emissions_price_so2)) %>%
  mutate(nonattain_emissions_price_so2 = ifelse(is.infinite(nonattain_emissions_price_so2), 0, nonattain_emissions_price_so2))

simulation_results <- simulation_results |>
  left_join(simulation_results_base) |>
  mutate(
    first_best_emissions_price_so2 = first_best_emissions_price_so2 - no_reg_emissions_price_so2,
    nonattain_emissions_price_so2 = nonattain_emissions_price_so2 - no_reg_emissions_price_so2,
  )

plot_sdf <- counties_sdf_simpler %>%
  left_join(simulation_results) %>%
  replace_na(list(nonattainment_1997 = FALSE)) %>%
  st_transform(5070) |>
  st_simplify(dTolerance = 1e3)

ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = (first_best_emissions_price_so2 + 1) / 1000, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "SO2\nFirst-Best\nEmissions\nPrice\n(thousand $)",
    colors = rev(c("#b2182b", "#f4a582", "#fddbc7", "#ffffff")),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    trans = "log10",
    labels = scales::label_dollar(),
    limits = c(NA, 100),
    breaks = c(.01, .1, 1, 10, 100)
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.96, 0.30)
  )

ggsave(
  paste0("output/first_best_prices_so2.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)
ggsave(
  paste0("output/figure-e5-d-first_best_prices_so2.eps"),
  width = 20,
  height = 10,
  dpi = 175
)

# first-best emission prices map nox

non_attainment <- fread("data/productivity_panel.csv") %>%
  filter(year >= 1986 & year <= 1997) %>%
  group_by(fips, industry) %>%
  mutate(total_nonattainment = sum(nonattainment, na.rm = T)) %>%
  ungroup() %>%
  mutate(nonattainment_1991 = nonattainment > 0 & year == 1991) %>%
  mutate(nonattainment_1997 = nonattainment > 0 & year == 1997) %>%
  group_by(fips) %>%
  mutate(nonattainment_1991 = max(nonattainment_1991)) %>%
  mutate(ever_nonattainment = max(nonattainment_1997) > 0) %>%
  ungroup() %>%
  filter(year == 1997) %>%
  mutate(fips = str_pad(fips, 5, "left", "0")) %>%
  distinct(fips, nonattainment_1997)

# Load the county data from the maps package
state_fips <- state.fips %>%
  distinct(fips, .keep_all = TRUE)

dup_names <- str_sub(state_fips$polyname, 1, str_locate(state_fips$polyname, ":")[, 1] - 1)

state_fips$polyname[!is.na(dup_names)] <- dup_names[!is.na(dup_names)]

simulation_results <- read.csv("data/counterfactual/welfare_opt-opt_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(nox_price_change = cf_emissions_price_nox / base_emissions_price_nox) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  filter(industry == "manufacturing") |>
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment) |>
  mutate(
    cf_emissions_price_nox = cf_emissions_price_nox * adjustment,
    base_emissions_price_nox = base_emissions_price_nox * adjustment
  ) |>
  select(fips, nox_price_change, nonattainment_1997, first_best_emissions_price_nox = cf_emissions_price_nox) %>%
  mutate(first_best_emissions_price_nox = ifelse(is.infinite(first_best_emissions_price_nox), 0, first_best_emissions_price_nox))

simulation_results_base <- read.csv("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  filter(industry == "manufacturing") |>
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment) |>
  mutate(
    cf_emissions_price_nox = cf_emissions_price_nox * adjustment,
    base_emissions_price_nox = base_emissions_price_nox * adjustment
  ) |>
  select(fips, nonattainment_1997, no_reg_emissions_price_nox = cf_emissions_price_nox, nonattain_emissions_price_nox = base_emissions_price_nox) %>%
  mutate(no_reg_emissions_price_nox = ifelse(is.infinite(no_reg_emissions_price_nox), 0, no_reg_emissions_price_nox)) %>%
  mutate(nonattain_emissions_price_nox = ifelse(is.infinite(nonattain_emissions_price_nox), 0, nonattain_emissions_price_nox))

simulation_results <- simulation_results |>
  left_join(simulation_results_base) |>
  mutate(
    first_best_emissions_price_nox = first_best_emissions_price_nox - no_reg_emissions_price_nox,
    nonattain_emissions_price_nox = nonattain_emissions_price_nox - no_reg_emissions_price_nox,
  )

plot_sdf <- counties_sdf_simpler %>%
  left_join(simulation_results) %>%
  replace_na(list(nonattainment_1997 = FALSE)) %>%
  st_transform(5070) |>
  st_simplify(dTolerance = 1e3)

ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = (first_best_emissions_price_nox + 1) / 1000, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "NOx\nFirst-Best\nEmissions\nPrice\n(thousand $)",
    colors = rev(c("#b2182b", "#f4a582", "#fddbc7", "#ffffff")),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    trans = "log10",
    labels = scales::label_dollar(),
    limits = c(NA, 100),
    breaks = c(.01, .1, 1, 10, 100)
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.96, 0.30)
  )
ggsave(
  paste0("output/first_best_prices_nox.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)

ggsave(
  paste0("output/figure-e5-b-first_best_prices_nox.eps"),
  width = 20,
  height = 10,
  dpi = 175
)

# ========================
# First-best prices
# ========================

# first-best emission prices map pm25

non_attainment <- fread("data/productivity_panel.csv") %>%
  filter(year >= 1986 & year <= 1997) %>%
  group_by(fips, industry) %>%
  mutate(total_nonattainment = sum(nonattainment, na.rm = T)) %>%
  ungroup() %>%
  mutate(nonattainment_1991 = nonattainment > 0 & year == 1991) %>%
  mutate(nonattainment_1997 = nonattainment > 0 & year == 1997) %>%
  group_by(fips) %>%
  mutate(nonattainment_1991 = max(nonattainment_1991)) %>%
  mutate(ever_nonattainment = max(nonattainment_1997) > 0) %>%
  ungroup() %>%
  filter(year == 1997) %>%
  mutate(fips = str_pad(fips, 5, "left", "0")) %>%
  distinct(fips, nonattainment_1997)



simulation_results <- read.csv("data/counterfactual/welfare_opt-opt_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(pm25_price_change = cf_emissions_price_pm25 / base_emissions_price_pm25) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  filter(industry == "manufacturing") |>
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment) |>
  mutate(
    cf_emissions_price_pm25 = cf_emissions_price_pm25 * adjustment,
    base_emissions_price_pm25 = base_emissions_price_pm25 * adjustment
  ) |>
  select(fips, pm25_price_change, nonattainment_1997, first_best_emissions_price_pm25 = cf_emissions_price_pm25) %>%
  mutate(first_best_emissions_price_pm25 = ifelse(is.infinite(first_best_emissions_price_pm25), 0, first_best_emissions_price_pm25))

simulation_results_base <- read.csv("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  filter(industry == "manufacturing") |>
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment) |>
  mutate(
    cf_emissions_price_pm25 = cf_emissions_price_pm25 * adjustment,
    base_emissions_price_pm25 = base_emissions_price_pm25 * adjustment
  ) |>
  select(fips, nonattainment_1997, no_reg_emissions_price_pm25 = cf_emissions_price_pm25, nonattain_emissions_price_pm25 = base_emissions_price_pm25) %>%
  mutate(no_reg_emissions_price_pm25 = ifelse(is.infinite(no_reg_emissions_price_pm25), 0, no_reg_emissions_price_pm25)) %>%
  mutate(nonattain_emissions_price_pm25 = ifelse(is.infinite(nonattain_emissions_price_pm25), 0, nonattain_emissions_price_pm25))

simulation_results <- simulation_results |>
  left_join(simulation_results_base) |>
  mutate(
    first_best_emissions_price_pm25 = first_best_emissions_price_pm25 - no_reg_emissions_price_pm25,
    nonattain_emissions_price_pm25 = nonattain_emissions_price_pm25 - no_reg_emissions_price_pm25,
  )

plot_sdf <- counties_sdf_simpler %>%
  left_join(simulation_results) %>%
  replace_na(list(nonattainment_1997 = FALSE)) %>%
  st_transform(5070) |>
  st_simplify(dTolerance = 1e3)

# summary pm25 errors
plot_sdf <- plot_sdf |>
  mutate(error = (nonattain_emissions_price_pm25 - first_best_emissions_price_pm25) / 1000)
summary(plot_sdf[plot_sdf$nonattainment_1997 == TRUE, "error"])
summary(plot_sdf[plot_sdf$nonattainment_1997 == FALSE, "error"])

# average PM2.5 first best price in attainment counties
plot_sdf |>
  filter(nonattainment_1997 == FALSE) |>
  summary()


ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = (first_best_emissions_price_pm25 + 1) / 1000, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "PM 2.5\nFirst-Best\nEmissions\nPrice\n(thousand $)",
    colors = rev(c("#b2182b", "#f4a582", "#fddbc7", "#ffffff")),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    trans = "log10",
    labels = scales::label_dollar(),
    limits = c(NA, 100),
    breaks = c(.01, .1, 1, 10, 100)
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.96, 0.30)
  )

ggsave(
  paste0("output/first_best_prices_pm25.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)
ggsave(
  paste0("output/figure-e5-c-first_best_prices_pm25.eps"),
  width = 20,
  height = 10,
  dpi = 175
)

# first-best emission prices map voc

# employment in the model is normalized to 1 so need to multiply by total employment because we are summing over all people harmed

simulation_results <- read.csv("data/counterfactual/welfare_opt-opt_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(voc_price_change = cf_emissions_price_voc / base_emissions_price_voc) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  filter(industry == "manufacturing") |>
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment) |>
  mutate(
    cf_emissions_price_voc = cf_emissions_price_voc * adjustment,
    base_emissions_price_voc = base_emissions_price_voc * adjustment
  ) |>
  select(fips, voc_price_change, nonattainment_1997, first_best_emissions_price_voc = cf_emissions_price_voc) %>%
  mutate(first_best_emissions_price_voc = ifelse(is.infinite(first_best_emissions_price_voc), 0, first_best_emissions_price_voc))

simulation_results_base <- read.csv("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  filter(industry == "manufacturing") |>
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment) |>
  mutate(
    cf_emissions_price_voc = cf_emissions_price_voc * adjustment,
    base_emissions_price_voc = base_emissions_price_voc * adjustment
  ) |>
  select(fips, nonattainment_1997, no_reg_emissions_price_voc = cf_emissions_price_voc, nonattain_emissions_price_voc = base_emissions_price_voc) %>%
  mutate(no_reg_emissions_price_voc = ifelse(is.infinite(no_reg_emissions_price_voc), 0, no_reg_emissions_price_voc)) %>%
  mutate(nonattain_emissions_price_voc = ifelse(is.infinite(nonattain_emissions_price_voc), 0, nonattain_emissions_price_voc))

simulation_results <- simulation_results |>
  left_join(simulation_results_base) |>
  mutate(
    first_best_emissions_price_voc = first_best_emissions_price_voc - no_reg_emissions_price_voc,
    nonattain_emissions_price_voc = nonattain_emissions_price_voc - no_reg_emissions_price_voc,
  )

plot_sdf <- counties_sdf_simpler %>%
  left_join(simulation_results) %>%
  replace_na(list(nonattainment_1997 = FALSE)) %>%
  st_transform(5070) |>
  st_simplify(dTolerance = 1e3)

ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = (first_best_emissions_price_voc + 1) / 1000, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "VOC\nFirst-Best\nEmissions\nPrice\n(thousand $)",
    colors = rev(c("#b2182b", "#f4a582", "#fddbc7", "#ffffff")),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    trans = "log10",
    labels = scales::label_dollar(),
    limits = c(NA, 100),
    breaks = c(.01, .1, 1, 10, 100)
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.96, 0.30)
  )

ggsave(
  paste0("output/first_best_prices_voc.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)
ggsave(
  paste0("output/figure-e5-e-first_best_prices_voc.eps"),
  width = 20,
  height = 10,
  dpi = 175
)

realloc_results <- read.csv("data/counterfactual/welfare_opt-opt_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(
    welfare = (exp(welfare - 1) - 1) * 100,
    amenities - amenities - 1,
    real_wage = real_wage - 1
  ) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  group_by(fips) %>%
  summarise(
    welfare_realloc = weighted.mean(welfare, base_labor),
    labor = sum(labor) / sum(base_labor) * 100
  ) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment)


simulation_results <- read.csv("data/counterfactual/welfare_opt-opt_prod-0.0_amen-on_abatement-on_transport-on_mobility-off_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(
    welfare = (exp(welfare - 1) - 1) * 100,
    amenities - amenities - 1,
    real_wage = real_wage - 1
  ) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  group_by(fips) %>%
  summarise(
    welfare_no_realloc = weighted.mean(welfare, base_labor),
    labor_no_realloc = sum(labor) / sum(base_labor) * 100
  ) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(realloc_results) %>%
  mutate(realloc_value = welfare_realloc - welfare_no_realloc)

plot_sdf <- counties_sdf_simpler %>%
  left_join(simulation_results) %>%
  replace_na(list(nonattainment_1997 = FALSE)) %>%
  st_transform(5070)

quantile(plot_sdf$realloc_value, seq(0, 1, .05), na.rm = T)
ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = realloc_value, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "Labor\nReallocation\nWelfare (%)",
    colors = c("#b2182b", "#f4a582", "#fddbc7", "#ffffff", "#92c5de", "#4393c3", "#2166ac"),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    limits = c(-3, 3),
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    values = scales::rescale(c(-3, -2, -1, 0, 1, 2, 3)),
    oob = scales::squish,
    breaks = c(-3, -1.5, 0, 1.5, 3)
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.96, 0.30)
  )

ggsave(
  paste0("output/welfare_map_realloc_value_opt.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)

ggsave(
  paste0("output/figure-e6-welfare_map_realloc_value_opt.eps"),
  width = 20,
  height = 10,
  dpi = 175
)
